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We discuss the thermodynamics of a non-relativistic gas of bosons with a local repulsive inter- 
action. In particular, we compute the temperature and density dependence of pressure, energy 
and entropy-density, superfluid and condensate- fraction, correlation length, specific heat, isother- 
mal and adiabatic compressibility and various sound velocities. The T ^ limit approaches the 
quantum phase transition, while the universal critical behavior of a classical second order phase 
transition in the 0(2) universality class determines the region around the critical temperature. Our 
non-perturbative flow equations based on exact functional renormalization cover all regions in the 
phase diagram. 



I. INTRODUCTION 

A gas of non-relativistic bosons with a repulsive point- 
like interaction is one of the simplest interacting statis- 
tical systems. Since the first experimental realization [l| 
of Bose-Einstein condensation (BEC) with ultracold 
gases of bosonic atoms, important experimental advances 
have been achieved, for reviews see [1, p|. Thermody- 
namic observables like the specific heat [5| or properties 
of the phase transition like the critical exponent v Q 
have been measured in harmonic traps. Still, the theo- 
retical description of these apparently simple systems is 
far from being complete. This is partly due to the failure 
of perturbation theory for several interesting regions in 
the phase diagram: the quantum critical behavior in the 
limit of zero temperature [7] or the critical behavior near 
the second order phase transition between the disordered 
and the BEC-phase. Also the generic infrared behavior 
in the BEC-phase, which is associated to a Goldstone 
boson, is non-perturbative 0. This issue becomes par- 
ticularly important in two- or one-dimensional systems. 

In this paper we employ non-perturbative flow equa- 
tions for an investigation of the thermodynamic quanti- 
ties for interacting bosons in three dimensions. In this 
approach the complications of the long-distance physics 
can be isolated and solved by a stepwise integration of 
the fluctuation effects. Already in a simple truncation we 
obtain a rather complete picture of the temperature and 
density dependence of many thermodynamic observables. 
In turn, if some of these observables can be measured 
precisely, and if the finite size effects due to the trap are 
properly taken into account, one may use such observ- 
ables as compressibility or sound velocity as a precise 
"thermometer" of the system. Simultaneous knowledge 
of two independent thermodynamic observables will be 
sufficient to determine the thermodynamic parameters - 
the temperature T and the density n. 

Furthermore, we study the dependence of the observ- 
ables on the interaction strength. If the interaction 
strength can be varied experimentally, for example by 
its dependence on a magnetic fields, many aspects of 
our computations can be tested by observation. In par- 



ticular, we have computed several non-analytic features 
near the phase transition. While the critical exponents 
and amplitude ratios in the immediate vicinity of the 
critical temperature Tc are given by universal laws, the 
amplitudes themselves as well as the approach to criti- 
cality and the location and size of the critical region in 
the phase diagram are non-universal. The dependence of 
these quantities on the interaction strength may become 
an important test for non-perturbative methods. 



II. METHOD 

The thermodynamic properties of a statistical system 
in thermal and chemical equilibrium, are described by 
the grand canonical partition function 



Z = Tre"'^'^"^"' 



(1) 



The trace operation includes a summation over all possi- 
ble states of the system including varying particle number 
N. We use (3 = 1/T with units where h = kB = I- The 
object flc is the thermodynamic potential of the grand 
canonical partition function. It has the differential 



dno = -SdT - Ndfi-pdV. 



(2) 



The partition function has a functional integral represen- 
tation 



D(l)e 



-s[<p] 



(3) 



where / Dip = JJ^^ -j / d(p{T, x) sums over all possible 
field configurations lp{t,x). The variable t is a periodic 
euclidean time coordinate in the range r S (0, [3) and x 
denotes the usual position coordinate. 

The microscopic action S[(f\ includes the Hamiltonian 
and the chemical potential. For nonrelativistic bosons 
with a pointlike interaction it is given by 



^*(a.-A-^V + -A(v5»n. (4) 
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We use energy units where 2M = 1, with M the mass of 
the atoms. Apart from the chemical potential and the ki- 
netic energy we include a repulsive pointlike interaction. 
The interaction strength A > determines the scatter- 
ing length a. For A = the functional integral in Eq. 
([3|) is Gaussian and can be done analytically, however for 
A > this is not possible any more. For small A, or more 
precisely small an^^^, one might try to make a perturba- 
tive analysis by expanding in A, but this often leads to 
infrared problems. Furthermore, we also want to cover 
the case of strong interactions. 

The method we use to determine the partition function 
in Eq. Q is of a different kind. We do not perform the 
functional integral in one step. Instead, we integrate out 
fluctuations with large momenta first and fluctuations 
with small momenta later. From a physical point of view 
it is not reasonable to include the effect of fluctuations up 
to an infinetly large momentum scale. A microscopic ac- 
tion as in Eq. ^ is an effective description of the physics 
at some momentum scale A (somewhat smaller than the 
inverse Bohr radius ciq^)- At larger momentum or for 
larger energies, the physics is modified. For example, one 
might have to include photons and electrons into the de- 
scription and at even larger energies quarks and gluons. 
All this is not visible anymore in the "microscopic ac- 
tion" ^ - the fluctuations with momenta < A^ have 
already been "integrated out". Starting from Eq. 
the computations of thermodynamics should only include 
fluctuations with momenta smaller than the "ultraviolet 
cutoff" A. 

On the other hand, any real system or any given 
physical observation has also an effective "infrared cut- 
off" scale fcph, such that only fluctuations in the range 
fcpjj < 9^ < A^ matter. It makes no sense to include fluc- 
tuations with a wavelength larger than the size of some 
macroscopic probe. Alternatively, fcph may also be set 
by an inverse characteristic wavelength of some measure- 
ment device. At the scale fcph we would like to dispose 
of an object similar to the action ([4]), with the differ- 
ence that now all fluctuations with > fcpj^ instead of 

> A^ are integrated out. Such an "average action" Q 
F/cpj^ averages out all structures with wavelength smaller 

then fcpj^^ and only retains the "macroscopic" information 
about physics at the momentum scale fcph. It will consti- 
tute a Landau- type theory for the macroscopic physics. 
Thermodynamic quantities can be easily derived from 
Ffeph, since no further fluctuations need to be taken into 
account. 

The basic challenge of a computation of thermodynam- 
ics from given microscopic laws is to build a bridge from 
the microscopic action S = Fa to the average action 
Fj^pi^ . This can be a difficult task, since Ffc^j^ niay be 
rather different from Fa. Typically, a large part of the 
precise microscopic information is lost in the averaging 
procedure. On the other hand, new phenomena due to 
collective effects may appear in Fph, which are not di- 
rectly visible in Fa. The basic idea of our approach is 
to make the infrared cutoff scale fc variable. The result- 



ing flowing actiori^ or scale dependent average action F^ 
interpolates continuously between the microscopic action 
for fc = A and the macroscopic action Tk^^^ . Lowering fc 
from A to fcph the fluctuation effects are included step- 
wise. In this sense it realizes the Wilsonian idea of renor- 
malization, even though the implementation of a sliding 
infrared cutoff leads to several important conceptual and 
technical differences as compared to the sliding ultravio- 
let cutoff investigated in the flrst approaches to functional 
renormalization [10]. 

The dependence of the flowing action F^^ on the in- 
frared cutoff scale fc obeys an exact functional differen- 
tial equation It can be solved approximately by 
a truncation of the most general functional form of the 
flowing action. Such truncations do not have to rely on 
the expansion in some small parameter as the interaction 
strength and can describe physical phenomena that are 
"non-perturbative" . 

The average potential Uk{p) obtains from Ffe[(^] by us- 
ing for the argument a homogeneous (and r-independent) 
fleld <p, with p = (p*(p. An evaluation of the potential at 
its minimum, C/min, yields directly the grand canonical 
partition function. 



-I3Q.G 



with 



Umin = Uk=oipo) with — C/fc=o(p) 



(5) 



(6) 



P=Pa 



The expectation value or order parameter po can be as- 
sociated with the condensate density. The average po- 
tential Uk{p) obeys the exact flow equation 



k^-^U.ip) 



{Pii + P22 + 2Rk)kdkRk 



jq{Pll+Rk){P22+Rk) 

We use the Matsubara formalism with 



p2 • 

A 2 



/ -TV, qo = 2nnT, n S N, / - / 7^' 



(7) 



,(8) 



In our truncation we approximate the momentum de- 
pendence of the inverse propagator for the radial and 
angular (Goldstone) mode by 



Pii(g) = At ^V(g^U' ^2pU" 
P22{q) = A(f + Vql + U' 
Pl2{q) = Sqo. 



(9) 



The flow equations for the quantities A, 5, and V can be 
found in ref. p^ . 

The effective potential U is related to the pressure by 



which has the differential 



dp = s dT + ndfi. 



(10) 



(11) 
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Here we use s = S/V for the entropy density and n = 
N/V for the particle density. The formal infinite volume 
limit corresponds to kp^ — 0. We can use our method to 
determine many thermodynamic quantities. Derivatives 
of U with respect to T and and /i are taken numerically 
by solving the flow equation for close enough values of T 
and fi. The numerical effort is reduced and the accuracy 
increased by using an additional flow equation for 



(12) 



explicit expressions for the flow equations can be found 
in (l^ . The approximation scheme we use in this pa- 
per is basically the same as the one used there. Since 
we use an infrared cutoff only for momenta but not for 
frequencies, the correct ultraviolet convergence for the 
sum of the Matsubara frequencies is not automatically 
obeyed for the flow equations. We have checked that all 
thermodynamic quantities discussed in this paper show 
a satisfactory convergence of the Matsubara sum, except 
for the pressure. In the flow equation for pk we set the 
frequency coefficients to their microscopic values S — 1, 
V = ior very large Matsubara frequencies |go| > ^uv 
For bosons with a pointlike repulsive interaction we 
found in [T3| that the scattering length is bounded by 
the ultraviolet scale a < 37r/(4A). This is an effect 
due to quantum ffuctuations similar to the "triviality 
bound" for the Higgs scalar in the standard model of 
elementary particle physics. For a given value of the di- 
mensionless combination an^/^ we cannot choose A/n^^^ 
larger then 37r/(4an^/'^). For our numerical calculations 
we use K/'n}/^ « 10. Other momentum scales are set by 
the temperature and the chemical potential. The low- 
est nonzero Matsubara frequency gives the momentum 
scale A|, = 2tiT. For a Bose gas with a — Q one has 
Tc/n2/3 « 6.625 such that Arjn^^^ « 6.45. The momen- 
tum scale associated to the chemical potential is A^ — /i. 
For small temperatures and scattering length one finds 



1/3 



fi « 87ran and thus A^/n 

We finally note that the thermodynamic relations for 
intensive quantities can only involve dimensionless ratios. 
We may set the unit of momentum by n^/^. The ther- 
modynamic variables are then T jr?!'^ and fxjv?!^ . The 
thermodynamic relations will depend on the strength of 
the repulsive interaction A or the scattering length a, and 
therefore on a "concentration" type parameter av}!^ . 



III. DENSITY, SUPERFLUID DENSITY, 
CONDENSATE AND CORRELATION LENGTH 

Let us start our discussion of the thermodynamic prop- 
erties with the density. In the grand canonical formalism 
it is obtained by taking the derivative of the thermody- 
namic potential with respect to 



dp I 
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FIG. 1: (Color online) Density in units of the scattering length 
no?' as a function of the (rescaled) chemical potential fjiO? . 
We choose for the temperatures Ta^ — 2 ■ 10~^ (solid curve), 
Ta^ = 4 ■ 10"'' (dashed-dotted curve) and Ta^ = 6 • 10"* 
(dashed curve). For all three curves we use aA = 0.1. 



We could compute the /i-derivative of p numerically by 
solving the ffow equation for U with neighboring values 
of /i. In we also describe another method which 
employs a flow equation directly for n. Since we often 
express dimensionful quantities in units of the interparti- 
cle distance n"^/'^, it is crucial to have an accurate value 
for the density n. Comparison of the numerical evalua- 
tion and the solution of a separate flow equation for n 
shows higher precision for the latter method and we will 
therefore employ the flow equation. We plot in Fig. [1] 
the density in units of the scattering length, na^ , as a 
function of the dimensionless combination fia? . 

For a comparison with experimentally accessible quan- 
tities we have to replace the interaction parameter A in 
the microscopic action by a scattering length a which 
is a macroscopic quantity. For this purpose we start the 
flow at the UV-scale Ayv with a given A, and then com- 
pute the scattering length in vacuum (T = n = 0) by 
following the flow to A; = [l^]- This is a standard pro- 
cedure in quantum field theory, where a "bare coupling" 
(A) is replaced by a renormalized coupling (a). For an 
investigation of the role of the strength of the interaction 
we may consider different values of the "concentration" 
c = av}/^ or of the product /xa^. While the concentration 
is easier to access for observation, it is also numerically 
more demanding since for every value of the parameters 
one has to tune /i in order to obtain the appropriate den- 
sity. For this reason we rather present results for three 
values of /za^, i. e. = 2.6 x 10"^ (case I), fxa? = 0.0040 
(case II) and fio? = 0.044 (case III). The prize for the nu- 
merical simplicity is a week temperature dependence of 
the concentration c — a'n}/'^ for the three different cases, 
as shown in Fig. [51 Here and in the following figures 
case I, which corresponds to an}/^ ~ 0.01, is represented 
by the little crosses, case II with an^^^ « 0.05 by the 
dots and case III with an^^^ ~ 0.1 by the stars. It is 
well known that the critical temperature depends on the 



concentration c 



an 



From our calculation we find 
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FIG. 2: (Color online) Concentration c = an^^^ as a function 
of temperature T/[rT?^^) for the three cases investigated in 
this paper. Case I corresponds to ain}^^ ~ 0.01 (crosses), case 
II to an^/^ « 0.05 (dots) and case III has an^^^ « 0.01 (stars). 




FIG. 3: (Color online) Superfluid fraction of the density ns/n 
as a function of the temperature T/n^^^ for the cases I, II, and 
III. For small T jr?^^ we also show the corresponding curves 
obtained in the Bogoliubov approximation (dashed lines). 



Tc/(n2/3) = 6.74 with c = 0.0083 at T = in case I, 
Tc/(n2/3) = 7.16 with c = 0.044 at T = in case II and 
finally T,J{n'^/'^) = 7.75 with c = 0.088 at T = in case 
III. 

This values can are obtained by following the super- 
fluid fraction of the density ns/n^ or equivalently the 
condensate part of the density nc/n as a function of 
temperature. For small temperatures T — > all of the 
density is superfluid, which is a consequence of Galilean 
symmetry. However, in contrast to the ideal gas, not all 
particles are in the condensate. For T = this conden- 
sate depletion is completely due to quantum fluctuations. 
With increasing temperature both the superfluid density 
and the condensate decrease and vanish eventually at 
the critical temperature T = Tc- That the melting of the 
condensate is continuous shows that the phase transition 
is of second order. We plot our results for the super- 
fluid fraction in Fig. [3] and for the condensate in Fig. |4l 
For small temperatures, we also show the corresponding 
result obtained in the framework of Bogoliubov theory 
[Tst (dashed lines). This approximation assumes a gas of 
non-interacting quasiparticles (phonons) with dispersion 
relation 

e{p) = V2Anp2+p4_ (14) 

It is is valid in the regime with small temperatures 
T Tc and small interaction strength an^/^ <^ i por 
a detailed discussion of Bogoliubov theory and the cal- 
culation of thermodynamic observables in this frame- 
work we refer to ref. [4]. Our curves for the superfluid 
fraction match the Bogoliubov result for temperatures 
T jr?^"^ < 1 in all three cases I, II, and III. For larger 
temperatures there are deviations as expected. For the 
condensate density, there is already notable a deviation 
at small temperatures for case III with av}!^ k, 0.1. This 
is also expected, since Bogoliubov theory gives only the 
first order contribution to the condensate depletion in a 
perturbative expansion for small an^l^ . For temper- 



nc/n 




FIG. 4: (Color online) Condensate fraction of the density 
nc/n as a function of the temperature T/n^^^ for the cases 
I, II, and III. For small T /n^^^ we also show the correspond- 
ing curves obtained in the Bogoliubov approximation (dashed 
lines) . 

atures slightly smaller than the critical temperature Tc 
one expects that the condensate density behaves like 

nc{T)=B^(^^^y (15) 

with (3 = 0.3485 the critical exponent of the three- 
dimensional XY-universality class [145 • Indeed, the con- 
densate density is given by nc — (t>o4'o where (j)Q is the 
expectation value of the boson field which serves as an 
order parameter in close analogy to e. g. the magnetiza- 
tion M in a ferromagnet. Eq. (I15|) is compatible with 
our findings, although our numerical resolution does not 
allow for a precise determination of the exponent (3. 

With our method we can also calculate the correlation 
length ^. For temperatures T < Tc one distinguishes 
between the Goldstone correlation length and the ra- 
dial correlation length S^j^. While the former is infinite, 
S,Q^ = 0, the latter is finite for T < Tc- It is also known 
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FIG. 5: (Color online) Correlation length ^ru^^'-^ for T < Tc 
and ^n^''^ for T > Tc as a function of the temperature T/n^^'^ 
for the cases I, II, and III. 



s/n 




FIG. 6: (Color online) Entropy per particle s/n as a function 
of the dimensionless temperature T jr?^'^ for the cases I, II, 
and III. For T/v?^^ < 5 we also plot the results obtained 
within the Bogoliubov approximation (dashed lines). 



as the "healing length" , given by 



Cr = 2Xpo = 2 



1 d^U 
A dp 



Po 



(16) 



and diverges only close to the phase transition. In the 
symmetric regime for T > Tc there is only one correlation 



length ^ 



^ ^ , which also diverges for T 



From the theory of critical phenomena one expects close 
to Tc the behavior 



Tc-T 

Tc 
T-Tc 

Tc 



for T <Tc 

for T < Tc. (17) 



The critical exponent v — 0.6716 [TJ] is again the one 
of the three-dimensional XY- or 0(2) universality class. 
We plot our result for the correlation length in units of 
the interparticle distance ^/jn^/-^ for T < Tc and ^n^/'^ 
for r > Tc as a function of the temperature T /n?/^ in 
Fig. El 



IV. ENTROPY DENSITY, ENERGY DENSITY, 
AND SPECIFIC HEAT 



The next thermodynamic quantity we investigate is the 
entropy density s and the entropy per particle s/n. We 
can obtain the entropy as 



dp I 



(18) 



We compute the temperature derivative by numerical dif- 
ferentiation, using flows with neighboring values of T and 
show the result in Fig. [6l For small temperatures our re- 
sult coincides with the entropy of free quasiparticles in 
the Bogoliubov approximation (dashed lines in Fig. [5]). 
As it should be, the entropy per particle increases with 



the temperature. For small temperatures, the slope of 
this increase is smaller for larger concentration c. We 
may consider a change of the volume 1^ by a change in 
the trap, for example increasing V by making the trap 
shallower. If no particles are lost, N = const., an increase 
of V corresponds to a decrease of n. Furthermore, if the 
change is adiabatic, the entropy and the entropy per par- 
ticle, s/ji, remains constant. For constant s/n the ratio 
T/r?!"^ = 7 is fixed, as given by Fig. [5] for given s/n. 
An adiabatic increase of the volume therefore induces a 
lowering of the temperature, T = ^r?!'^ ~ 'V~'^l'^ . This 
allows for a continuous reversible variation of the temper- 
ature by means of a variation of the trap parameters. By 
an adiabatic increase of a (at fixed iV) one can increase 
the ratio T/n^/^. This can be realized by a variation of 
a magnetic field, which may therefore be used to explore 
the phase transition and the region of T /r?!"^ around the 
critical temperature. 

From the entropy density s we infer the specific heat 
per particle, 



T as 
n 9r 



(19) 



as the temperature derivative of the entropy density at 
constant particle density. Using the Jacobian, we can 
write 



ds I 



d{s,n) d{s,n)d{T,^) 



a(r, n) d{T, n) 5(T, n) ' 
For the specific heat this gives 




5s I dn . 



dn I 



(20) 



(21) 



Our result for the specific heat per particle is shown for 
different scattering lengths in Fig. [71 While this quantity 
is positive in the whole range of investigated tempera- 
tures, it is interesting to observe the cusp at the critical 



6 




FIG. 7: (Color online) Specific heat per particle c„ as a func- 
tion of the dimensionless temperature T/rt^^^ . The dashed 
lines show the Bogoliubov result for c„ which coincides with 
our findings for small temperature. However, the characteris- 
tic cusp behavior cannot be seen in a mean-field theory. 



p/(n^/3) 




FIG. 8: (Color online) Pressure in units of the density p/rC'^^ 
as a function of temperature T/v?^^ for the cases I (crosses), 
II (dots), and III (stars). We also show the curves ob- 
tained in the Bogoliubov approximation for small tempera- 
tures (dashed lines). 



temperature Tc which is characteristic for a second or- 
der phase transition. This behavior cannot be seen in a 
mean-field approximation, where fluctuations are taken 
into account only to second order in the fields. Only for 
small temperatures, our curve is close to the Bogoliubov 
approximation, shown by the dashed lines in Fig. [T] 

In fact, close to Tc the specific heat is expected to 
behave like 



hi - 6, 



Tc 
T-Tc 
Tc 



for T <Tc 



for T>Tc, (22) 



with the universal critical exponent a of the 3- 
dimensional XY universality class, a = —0.0146(8) p^ . 
The critical region, where the law Cy ~ \T — T^^" holds, 
may be quite small. Our numerical differentiation proce- 
dure cannot resolve the details of the cusp. 

In the grand canonical formalism, the energy density e 
is obtained as 



e = —p + Ts + fin. 



(23) 



(The total energy density density in the gravitational 
context reads nMc^ + e.) We plot p/{'rv'/'^) as a func- 
tion of temperature in Fig. [8] and the energy density 
t/{rr'/^) is plotted in Fig. O We have normalized the 
pressure such that it vanishes for T = /i = 0. Technically 
we subtract from the flow equation of the pressure the 
corresponding expression in the limit T = = 0. This 
procedure has to be handled with care and leads to an 
uncertainty in the offset of the pressure, i. e. the part 
that is independent oiT/m?/^ and [ijr?!'^. 

For zero temperature, the pressure is completely due 
to the repulsive interaction between the particles. For 
nonzero temperature, the pressure is increased by the 
thermal kinetic energy, of course. For the energy and the 




2/3 



Tin 



FIG. 9: (Color online) Energy per particle e/'n?^^ as a function 
of temperature T/'n?^^ for the cases I (crosses), II (dots), and 
III (stars). We also show the curves obtained in the Bogoli- 
ubov approximation for small temperatures (dashed lines). 



pressure we find some deviations from the Bogoliubov 
result already for small temperatures in cases II and III. 
These deviations may be partly due to the uncertainty 
in the normalization process described above. For weak 
interactions av}/^ ~ 0.01 as in case I, the Bogoliubov 
prediction coincides with our result. 



V. COMPRESSIBILITY 

The isothermal compressibility is defined as the rela- 
tive volume change at fixed temperature T and particle 
number N when some pressure is applied 



Kt — —T7 



1^1 



1 dn I 
n dp 



(24) 
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Very similar, the adiabatic compressibility is 



]_dV, 



1 dn I 
n dp 



(25) 



where now the entropy S and the particle number N are 
fixed. Let us first concentrate on the isothermal com- 
pressibility kt- To evaluate it in the grand canonical for- 
malism, we have to change variables to T and fi. With 
dp/diJ,\ rp — n and dp/dn\ — ndfi/dn\ one obtains 



1 dn I 



(26) 



This expression can be directly evaluated in our formal- 
ism by numerical differentiation with respect to /i. 

The approach to the adiabatic compressibility is simi- 
lar. Using again the Jacobian we have 



We need therefore 



1 dn, 1 d{n,s/n) 

nQp\s/n nd{p,s/n) 
ld{n,s/n) di^i,T) 
n d{pL,T) d{p,s/n) 



d{n, s/n) 
and also 



1 f dn . ds . dn . ds , 
n \d^^T^\^, ~ df^t^'d^J'^ 



(27) 



(28) 



dip, s/n) 


I / dp , ds : 


dp 


s dn 1 


diti,T) 


n\dfl^'^df^i' 


d^ 






dp , ds , 


dp 1 


s dn. \ 








n d^ ''^ J 



(29) 



^ 9s I ^s dn : 



n- 



dn 

^dJl^T 



In the last equations we used the Maxwell identity 
^1 = -i^L. Combining this we find 

dT Ifj, dp, IT 



[ dn I ds I f dn I \ j 

( 



n2ds_ I 



2 STi Q'J~' 



I „2dn\ 



(30) 



(9^p/9T^)|^ we need to evaluate a see- 



Since ds/dT\^ 

ond derivative numerically. We plot the isothermal and 
the adiabatic compressibility in Figs. [TUlandfTT] 

For the isothermal compressibility the temperature de- 
pendence is qualitatively different than in Bogoliubov 
theory already for small temperatures, while there seem 
to be only quantitative differences for the adiabatic com- 
pressibility. The perturbative calculation of the com- 
pressibility is difficult since it is diverging in the non- 
interacting limit an^/'^ 0. 




FIG. 10: (Color online) Isothermal compressibility kt n^^^ as 
a function of temperature T jr?^'^ for the cases I (crosses), II 
(dots), and III (stars). We also show the Bogoliubov result 
for small temperatures (dashed lines). 
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FIG. 11: (Color online) Adiabatic compressibility ks rc"^^ as 
a function of temperature T/n^^^ for the cases I (crosses), II 
(dots), and III (stars). We also show the Bogoliubov result 
for small temperatures (dashed lines). 



VI. ISOTHERMAL AND ADIABATIC SOUND 
VELOCITY 

The sound velocity of a normal fluid under isothermal 
conditions, i. e. for constant temperature T is given by 



1 dp I 



(31) 



We can obtain this directly from the isothermal com- 
pressibility 



{nnrY 



(32) 



as follows from Eq. (|24p . We plot our result for v'^ in Fig 
[T2I recalling our units 2M 



1 such that v\ stands for 
2Mv\. This plot also covers the superfluid phase where 
the physical meaning of v\ is partly lost. This comes 
since the sound propagation there has to be described by 
more complicated two-fluid hydrodynamics. In addition 
to the normal gas there is now also a superfluid fraction 



8 




FIG. 12: (Color online) Isothermal velocity ol sound as ap- 
propriate for single fluid v%lri'^^ = l/(«Tn^''"^) as a func- 
tion of the dimensionless temperature T jri'^^ for the cases I 
(crosses), 11 (dots), and 111 (stars). We also show the Bogoli- 
ubov result for small temperatures (dashed lines). 



2/3 




2/3 



Tjn 



FIG. 13: (Color online) Adiabatic velocity of sound as ap- 
propriate for single fluid v%/n^^^ = 1/ {us n^^^) as a func- 
tion of the dimensionless temperature T /v?^^ for the cases I 
(crosses), 11 (dots), and 111 (stars). We also show the Bogoli- 
ubov result for small temperatures (dashed lines). 



allowing for an additional oscillation mode. We will de- 
scribe the consequences of this in the next section. 

For most applications the adiabatic sound velocity 
is more important then the isothermal sound velocity. 
Keeping the entropy per particle fixed, we obtain 



and therefore 



^2 ^ J_^| 



Mvl = {nnsY 



(33) 



(34) 



Our numerical result is plotted in Fig. [131 Again the 
plot covers both the superfluid and the normal part, but 
only in the normal phase the object w| has its physical 
meaning as a sound velocity. 



VII. FIRST AND SECOND VELOCITY OF 
SOUND 

For temperatures < T < Tc there are two compo- 
nents of the gas: the superfluid and the normal part, ft 
was shown by Landau [1^ that this leads to two-fluid 
hydrodynamics with two distinct velocities of sound Ci/2 
corresponding to different kinds of excitations. 

The main reason for the existence of two sound veloci- 
ties is that the entropy flow is carried only be the normal 
component while the particle flow (or equivalently mass- 
flow) is carried by both the normal and the superfluid 
part. The continuity equation for the conserved particle 
number reads 



5(71 + V • j = 0, 



(35) 



where j ~ tlnVn + nsvs is the (complete) particle 
number current and -CV, vs are the velocities of the 
normal (tt-at) and superfluid (ng) parts of the density, 
n = Um + ns- The conservation equation for the entropy 
reads 



dts + V ■ [svn) = 0. 



(36) 



We work in linear order in an expansion in the velocities 
vn and vs- To close the set of hydrodynamic equations 
for small i/at, vs we need the equations for momentum 
conservation 

Mdt]+Vp^O, (37) 

and for the change in the superfluid velocity 

Mdtvs + V^l = Q. (38) 

The last equation guarantees that the superfluid flow re- 
mains irrotational, V x z/g = 0. 

From the combination of Eq. ((35)) and pT]) one obtains 

(39) 



Md'^n = Ap. 



To linear order in vs and vn one infers from the combi- 
nation of Eq. dSSl) and §^ 



nsV ■ {m - Vs) 



-dt{s/n). 



(40) 



We recover s/n = const, for = as appropriate for 
the disordered phase. Similarly, the combination of Eq. 
([37)1 and ^ gives 

MnNdt{vN — vs) ~ nV^ — \7p 

= -sVT. (41) 

The last equation uses the relation 

Vp = sVT + nV/i (42) 
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which follows directly from the differential of p, Eq. ([TT|) . 
Combining now Eqs. and (|^T|) yields the analogue 

ofEq. (pjl. 



2/3 



Mdf{s/n) = 4 — AT. 



(43) 



One next makes an ansatz for the thermodynamic vari- 
ables in the form 

p = Po + Sp, T = To + 5T 
n = no + Sn, s/n =^ sq/hq + d{s/n), (44) 

where po, Tq, no and sq are constant in space and time 
whereas 6p, 6T, Sn, and 5{s/n) are small and vary like 
sin[p{x—ct)] . We use ST and Sn as independent variables, 
with 



9p 

9ri 



0(s/7l) = — — 01 - 



Sn, 

d{s/n) 
dn 



, Sn, 



(45) 




FIG. 14: (Color online) First velocity of sound cl/n^^^ as 
a function of the dimensionless temperature T jr?^^ for the 
cases I (crosses), II (dots), and III (stars). We also show the 
prediction from Bogoliubov theory for T — > (short solid 
lines) . 



in order to obtain from Eqs. ([391) and (|43)l the wave 
equation 



2/3 



Mc 



2 d{s/n 



_ s ns 

dT I n n'-^riN 

dp I 

dT \n 



on IT 
on It 




= 0. 



(46) 



As a condition for possible sound velocities c one obtains 



dp 



dn 



+ 



s^nsT 



s/n 



n^npfC-u 



+ 



s^nsT dp 



n'^n^Cv dn 



(47) 




T/n 



2/3 



This relation uses 



as well as 



T 



d{s/n) . 
dT I 



dis/n) 


CudT. 


dp , dp , 




dn 


T T dp^'^ 


dn''^ dn'^/^ 




The latter relation follows from 




dp 


dp 1 


dp: dT. 




dn 


s/n = -q:^\t ^ 






together with 








dT: 


d{T,s/n) _ 


d{T,s/n) d{T,n) 




d(n, s/n) 


d(T,n) d{s/ 


n, n) 



(48) 



(49) 



(50) 



d{s/n), T 
dn ^ Cv 



(51) 



With these ingredients one can now solve Eq. P7|) for 
the first and second velocity of sound. The numerical 



FIG. 15: (Color online) Second velocity of sound c^/ri^^'"'' as 
a function of the dimensionless temperature T /r?^^ for the 
cases I (crosses), II (dots), and III (stars). We also show the 
prediction from Bogoliubov theory for T — * (short solid 
lines). For T /v?^^ < 1 our numerical determination becomes 
unreliable, since C2 is dominated by a ratio of terms that 
vanish for T ^ 0. 



results as a function of temperature are shown in Fig. 
[141 and [151 We also show there the prediction from 
Bogoliubov theory for T ^ (short solid lines). For 
cf the agreement with our findings is rather good, al- 
though there are some deviations for strong interactions 
as in case III. For C2 our numerical determination be- 
comes unreliable for T/n^^'^ < 1 since C2 is dominated 
by the term {s'^nsT) / {n^nNCy) in Eq. (|47l) . In the limit 
T ^ the quantities s, njv, and Cy also go to zero so 
that the numerical value for C2 is sensitive to the precise 
way how this limit is approached. 



10 



We observe that Eq. p?)) can be written as 



2\2 



{n ~ ns)cv-n? 



0, 



(52) 



with the single fluid isothermal and adiabatic sound ve- 
locities vt and vs given by Eqs. ((3T|) and ((33|) . This 
shows that c coincides with ws in the disordered phase 
where ns = 0- An intuitive form of the wave equation 
can be written as 



d^6n = vlA6n+ (vl -v^)ASf, 
d^Sf = {vl-vl + v^)ASf + vlASn, 



with 



and 



V 



s^nsT 



ii?{n - ns)cy ' 
Td{s/n). 



dn 



dT. 
dn 



(53) 



(54) 



(55) 



For fluctuations of 5n and 5T only vt, vs and v mat- 
ter. In the limit T ^ one observes Vg such that 
the fluctuations Sn are governed by the isothermal sound 
velocity vt- On the other hand, the the velocity v char- 
acterizes the dynamics of a linear combination of ST and 
dn. 



perature, density and the interaction strength of vari- 
ous thermodynamic observables: entropy, pressure, en- 
ergy, specific heat, isothermal and adiabatic compress- 
ibility, isothermal and adiabatic sound velocity as well as 
the first and second velocity of sound in the superfiuid 
phase. Non-analytic features at the critical temperature 
are clearly visible even though the limits of numerical 
differentiation do not allow a very fine resolution. The 
truncation of the functional flow remains rather simple, 
but we do not expect qualitative changes from the use 
of an extended truncation. The present truncation is 
already sufficient for reproducing correctly all expected 
qualitative features, both near the critical temperature 
of the phase transition between the superfiuid and dis- 
ordered phase and the quantum phase transition in the 
zero temperature limit. The long distance physics is un- 
der control and no infrared problems occur. Quantitative 
improvements may be achieved by an extension of the 
truncation and by increased numerical precision near the 
critical temperature. 

The computation of the thermodynamic response func- 
tions enters the hydrodynamic equations. They can be 
used for an investigation of the motions of atoms in a 
trap. Precise knowledge of the thermodynamics may 
allow for precision studies of this motion under the in- 
fluence of time varying trap geometry or time varying 
coupling strength. 



VIII. CONCLUSIONS 
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